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We present a model on the simulation of the measurement-device independent quantum key distribution 
(MDI-QKD) with phase randomized general sources. It can be used to predict experimental observations of 
a MDI-QKD with linear channel loss, simulating corresponding values for the gains, the error rates in 
different basis, and also the final key rates. Our model can be applicable to the MDI-QKDs with arbitrary 
probabilistic mixture of different photon states or using any coding schemes. Therefore, it is useful in 
characterizing and evaluating the performance of the MDI-QKD protocol, making it a valuable tool in 
studying the quantum key distributions. 

There has been a long history between the attacks and the anti-attacks in the development of quantum key 
distributions (QKD) since the idea of BB84 (Bennett-Brassard 1984'"^) protocol was put forward, due to the 
conflictions between the "in-principle" unconditional security and realistic implementations. Till today, 
there have been many different proposals for the secure QKD with realistic setups, such as the decoy-state 
method*""^' which can rescue the QKD with imperfect single-photon sources'"* while the device-independent 
quantum key distribution"" and the recently proposed measurement- device independent quantum key distri- 
bution (MDI-QKD)'"'"'' can further relieve the QKD even when the detectors are controlled by the eavesdrop- 
pgj.40,41 ]y[Q5^ interestingly, the MDI-QKD is not only immune to any detector attacks, but also able to generate a 
significant key rate with existing technologies. Moreover, its security can still be maintained with imperfect single- 
photon sources'"'"''' ""'^ 

In developing practical QKDs, one important question is how to evaluate the performance of a proposal before 
really implementing it, since it is not realistic to experimentally test everything. Therefore, it is an important job to 
make a theoretical study and numerical simulation to predict the experimental results. In principle, it allows to use 
different kinds of sources in a decoy state MDI-QKD''^''". Before experimentally testing all of them, one can choose 
to give a theoretical comparison with a reasonable model. In the traditional decoy state methods* the simu- 
lation model with linear channel loss is relatively simple. For example, given a source state, it is easy to calculate 
the state after transmission, and further estimate the gains and error rates possibly being observed in an experi- 
ment. However, for MDI-QKDs, it is not a simple job except for the special case of using weak coherent states, 
because both sides send out pulses and the successful events and errors are defined differently. So far, there have 
been proposals with different sources, e.g., the heralded single-photon source (HSPS) etc^"^'*^. And it has been 
shown that such a source can promise a longer secure distance than the weak coherent state. Nevertheless, it is 
unknown whether there are other sources which can present even better performance. Therefore, a general model 
on simulating the performance of arbitrary source states will be highly desirable. Here in this manuscript we solve 
the problem. 

For simplicity, we assume a linear lossy channel in our model. Note that the security does not depend on the 
condition of linear loss at all. We only use this model to predict: what values the gains and error rates would 
possibly be observed if one did the experiment in the normal case when there is no eavesdropper. Given these 
values, one can then calculate the low bound of the yield and the upper bound of the phase flip-error rates for 
single-photon pairs. The major goal here is to simulate the values of gains and error rates of different states in 
normal situations. Of course, they can be replaced with the observed values in real implementations. 

The paper is arranged as follows: In Sec. II we present the general model for the gains and error rates in a MDI- 
QKD, describing the detailed calculation processes. In Sec. Ill we proceed corresponding numerical simulations. 
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comparing the different behaviors of MDI-QKDs when using differ- 
ent source states. Finally, discussions and summaries are given out in 
Sec. IV. 

The general model on MDI-QKD 

A. Setups and definitions. Consider the schematic setup in Fig. P", 
there are three parties, the users- Alice and Bob, and the un-trusted 
third party (UTP)-Charlie. Alice and Bob send their polarized 
photon pulses to the UTP who will take collective measurement on 
the pulse-pairs. The collective measurement results at the UTP 
determine the successful events. They are two-fold click of 
detectors (1,4), (2,3), (1,2) or (3,4). The gain of any (two-pulse) 
source is determined by the number of successful events from the 
source. There are 4 detectors at the UTP, we assume each of them has 
the same dark count rate d, and the same detection efficiency C. In 
such a case, we can simplify our model by attributing the limited 
detection efficiency to the channel loss. Say, if the actual channel 
transmittance from Alice to Charlie is jji, we shall assume perfect 
detection efficiency for Charlie's detectors with channel transmit- 
tance oir\i^. Each detector will detect one of the 4 different modes. 



say a'^, Uy, bf^, by in creation operator. For simplicity, we denote 
them by c], i.e., c| = flj^, C2 = a\,, c| = &]^, c\ = by. In such a way, 
detector D, corresponds to mode i exactly. To calculate the gains 
that would-be observed for different source states in the linear 
lossy channel, we need to model the probabilities of different 
successful events conditional on different states. Let's first 
postulate some definitions before further study. 

Definition 1: event {i,j). We define event as the event that both 
detector and detector j click while other detectors do not click. 
Obviously, each j must be from numbers {1, 2, 3, 4} and i j. 
For simplicity, we request < j throughout this paper. 

Definition 2: Output states and conditional probabilities of each 
events: notations Pg^{. the output state of the beam-splitter. |/,-, /y) = 
I lilj): the beam-splitter's specific output state of /, photon in mode i, Ij 
photon in mode j, and no photon in any other mode. Explicitly, 



1'.';) 



' |0). P{i}\li, Ij) and P{ij\Pout)- the probability that 



event (i, j) happens conditional on that the beam-splitter's output 
state is |/,Z^) and Po„(, respectively. Hereafter, we omit the comma 
between Z, and Ip i.e., we use \ ljlj) for /,), and P{ij\lilj) for P{ij\li, Ij). 

Definition 3: Events' probability conditional on the beam-splitter's 
input state: p^!\ki,k2) =p'^j\kik2).We denote p^'^(fci,fc2) as the prob- 
ability of event (;, j) conditional on that there are ki photons of 
polarization a for mode a and /c2 photons of polarization P for mode 
b as the input state of the beam-splitter. Hereafter, we omit the 



PBS 




PBS 



Figure 1 | A schematic of the experimental setup for the collective 
measurements at the UTP. BS: beam-splitter; PBS: polarization beam- 
splitter; Dl - D4: single-photon detector; 1, 2: input port for photons. 
Incident and output lights at the left side of BS are mode a, lights at the 
right side of BS are mode b. 



comma between ki and ^2. ix or /? indicate the photon polarization. 
Explicitly, a or /J can he H,V, +, — for polarizations of horizontal, 
vertical, n/4 and 37:/4, respectively. To indicate the corresponding 
polarization state, we simply put each of these symbols inside a ket. 

Definition 4: Events' probability conditional on the two-pulse state 
of Alice and Bob's source: c{^{Pa®Pb)-^^ is the probability that event 
(i, f) happens conditional on that Alice sends out photon-number 
state Pa with polarization a and Bob sends out photon number state 

with polarization /J. Sometimes we simply use g^''* for simplicity. 

B. Elementary formulas and outline for the model. Given the defi- 
nitions above, we now formulate various conditional probabilities. 
We start with the probability of event {i,f) conditional on the output 
state \ lilj). 

{{l-df, if /,>0, lj>Q 
d{\-df, if l,-lj = 0 and l, + lj>0 (1) 
d^l-df,i{ l, = lj^O 

Here the detection efficiency does not appear because we put shall 
this into the channel loss and hence we assume perfect detection 
efficiency. The factor (1 — dy comes from the fact that we request 
detectors other than /, j not to click. Also, the probability for event 
{i,j) is 0 if any mode other than i,j is not vacuum. Given these, we can 
now calculate probability distribution of the various two fold events 
given arbitrary input states of the beam-splitter. Therefore, for any 
output state of the beam-splitter Po„,, the probability that event (i, j) 
happens is 

p{ij\Pout)-J2piWih)mPout\hij) (2) 

Based on this important formula, we can calculate the probability of 
event (/, j) for any input state by this formula. For the purpose, we 
only need to formulate Po„(. Therefore, given the source state of the 
two pulses Pa ® Ps-, we can use the following procedure to calculate 
the probability of event (i,j),pijipA ® Pb)' 

i) Using the linear channel loss model to calculate the two-pulse 
state when arriving at the beam-splitter. Explicitly, if the chan- 
nel transmittance is )/, any state j«) («| is changed into 

i«)(«HEcy(i-'i)""V)(fci- (3) 



ii) Using the transformation: a'^j y - 

— ^ (^alj V • to calculate the output state of the beam- 

splitter, pout- 

ill) Using Eq.(2) to calculate the probability of event (;', j). 
According to the protocol, we shall only be interested in the 
probabilities of successful events, (1, 2), (3, 4), (1, 4) and (2, 3). 
Below we will describe the detailed calculation processes in Z 
basis and X basis individually. 

In Z basis, all successful events correspond to correct bit values 
when Alice and Bob send out orthogonal polarizations, and they 
correspond to wrong bit values when Alice and Bob send out the 
same polarizations. The observed gain in Z basis for photon-number 
state Pa ® Pb is. 



^i^Ps = J J2 UfipA®PB)+q^"ipA®PB) + 



{i,j)eSuc 

^"{pA®PB)+qY{pA®PB: 



(4) 



and the set Sue = {(1, 2), (3, 4), (1, 4), (2, 3)). Here, as defined in 
Definition 4, q^f{pA®PB) represents the probabOity of event (;, j) 
whenever Alice sends her photon number state pA with polarization 



SCIENTIFIC REPORTS | 4:4612 | DOI: 1 0. 1 038/srep0461 2 



2 



a and Bob sends his photon number state pg with polarization For 
simplicity, we shall omit ® pg in brackets or in subscripts if there 
is no confusion. Meantime, the successful events caused by the same 
polarizations will be counted as wrong bits. These will contribute to 
the bit-flip rate by: 



4SZ 



(5) 



In X basis, we should be careful that the situation is different from in 
Z basis, since now the successful events correspond to correct bits 
include two parts: 1) Alice and Bob send out the same polarizations 

(+ + or ), and Charlie detects <1>* ((1,2) or (3,4) events happen); 

2) Alice and Bob send out orthogonal polarizations (H — or — h), 
and Charlie detects 'P" ((1,4) or (2,3) events happen). And the left 
successful events belong to wrong bits. Therefore, we have 



E 

{i,j)eSuc 



,+ + . 

V 



(6) 



and 



,^_ E(,-.;)E(14).(23) ] +E(i,j)E(12).(34) [^.f + lij 

4SX (7) 

Moreover, there are alignment errors which will cause a fraction (£^) 
of states to be flipped. We then modify the error rate in different bases 
by 

E^=Ed-(l-2E^)+E^ (8) 



and 



E^ = Ed-(l-2¥)+E^ 



(9) 



Note that in the above two formulas above, we have considered this 
fact: before taking the alignment error into consideration, the suc- 
cessful events can be classified into two classes: one class has no error 
and the other class has an error rate of 50%, they are totally random 
bits. The second class takes a fraction of 2E^ (or 2E'*") among all 
successful events. Alignment error does not change the error rate 
of the second class of events, since they are random bits only. 

Given these, we can simulate the final key rate. In a model of 
numerical simulation, our goal is to deduce the probably would-be 
value for S^, and E^, E^ in experiments. One can then calculate the 
yield of the single-photon pairs, Sn, the bit-flip rates in Zbasis andX 
basis, and hence the final key rate. Now everything is reduced to 
calculate all pjj'* above. 

C. Conditional probabilities for beam-splitter's incident state of 
fei photons in mode a and kj photons in mode b. We consider the 
case that there are incident photons in mode a and ^2 incident 
photons in mode b of the beam-splitter. Each incident pulse of the 
beam-splitter has its own polarization and is indicated by a subscript. 
In general, we consider the state 

(10) 

We shall consider the conditional probabilities for various successful 
events, i.e. p^j\kik2). Since we only consider the incident state of ki 
photons in mode a and ^2 photons in mode b, we shall simply usepjj'* 
forp^^(/ciA:2) in what follows. 

i) in Z basis 

First, we consider the following two-mode state 

1 



\h)H\k2)v = 



11 



as the input state of the beam-splitter. After BS, the output state | xj/) is 



Therefore 



ki+ki 



a'y-bU |0) (12) 



(/i/2|p„„,|/i/2) = (l/2)'^'+'^^<5,,,.<5,,t (13) 
According to Eq.(2), the conditional probability for event (1,2) is 



p'^,^=p{u\p„j=Y,pm,h){i/2f^'^5,,i,5,j, 



(14) 



= (l/2)'''+'^^P(l2|/ci/c2) 
Similarly, we have 



o^y={i/2f+'^P(M\hk2) 



\ki +k2 , 



pf/ = (l/2)^'-^^P(14|fc,fc2) 



(15) 



p^^ = {l/2r'^'^P{23\k2h) 



Note that here P{ij\k^k„) is just P{ij\li = k^, Ij = k„) when li = ki as 
defined by our Definition 2 in previous section. For example, 
P{2?)\k2ki) is actually P(23|/2 = kj, I3 = ki). Similarly, if the beam- 
splitter's input state is \ki)y\k2)H> i-e- vertical photons in mode a 
and kz horizontal photons in mode b, we have 

prf = (l/2)*' + '^P(12|/c2fcO 

P3^« = (l/2)*^' + '^^P(34|/c2fci) 

pri' = (l/2)*' + '^P(14|fc2fci) 

p™ = {l/2f + '^P{23\kik2) 

Next we consider the following two-mode state 

1 



(16) 



\h)H\k2)H- 



al''bl''\0} 



(17) 



as the input state of the beam-splitter. After the beam-splitter, it 
changes into: 

hA>^(^)'''^(4+^I.)'(«I,-;^l.)"|o) (18) 

We have the following uniform formula for probabilities of any 
successful events: 



k,\k2\ 

{ ki + h) 

k,\k2\ 



(l/2)'"+*'P(y|A:i -1-^2,0); fori = l,7 = 2; ori = 3,j = 4 



(l/2)'" + ''^P{ij\0,ki+k2); for! = l,; = 4; or! = 2,; = 3 

(19) 

Similarly, when the beam-splitter's input pulses are both vertical, we can 
find the value forpj^. 

ii) X basis 

We first consider the beam-splitter's input state of jA:i)+|fc2)-> i-e-i 
there are k^ photon with n/4 polarization in mode a and kz photons 

with 37t/4 polarization in mode b. Note that | + ) = {\E[) ±\V)). 

v2 

The output state of the beam-splitter is 



1 



'al-al-bl + blf\0) 



i+ai+bi+biy^ 



(20) 
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We have 



where 



Ai = min{lj,ki}, A2 = — min{lj,k2} 



+ k2 

(21) 



(22) 



and min{li, fci(fc2)} is the smaller one of /, and ^1(^2). Thus we can 
calculate the conditional probabilities by 

l, = 0 



Hence 



A, 



P(ij\l,M+k2-l,) 
for / = l,j = 2 and / = 3, j = 4; and 

P(!)l/„fci+fc2-/,) 



^ V/,!(fci+fc2-;,)!a.c^r(-i; 



(23) 



£ J2 VkKh +k2-ii)icid, 



for ! = 1, j = 4 and i = 2,j = 3. Besides, it is easy to show 



(24) 



(25) 



If the polarization of incident pulses of the beam-splitter are both 
7i/4, then the output state is 



2i"+i'^^/kjk^. 



al + al + bl + bl 



(26) 



al + a\-hl-b\) » 



We find 



'-' ^' (,=0 s = A2 

P(ii\l,M+k2-li) 
for ! = l,j = 2 and i = 3,j = 4; and 



(27) 



^'^^"4^^+^^Jtilfc2! ^ 

P(!)l/„fci+fc2-;,) 



Ai 



^ v//,!(/c,+fc2-/,-)!a.c^r^(-i)'" 



for ! = 1, j = 4 and i = 2,j = 3. Also, we have 



(28) 



(29) 



D. ProbabUities of events conditional on source states. In the above 
subsection, we have formulated the probabilities of various events 
conditional on a pure input state jfci)|fc2). In fact, the results can be 
easily extended to the more general case when the beam-splitter's 
input state is a mixed state. Say, 



Y^fkl\kl){h]^[J2f,2\k2){k2 



(30) 



Suppose the polarizations of mode a, b are v., /!, respectively. We then 
have 



t^f=Y.f^Apf(hk2) 



(31) 



where pJJ^ (fci A:2 ) is the same as defined in the previous subsection, for 
all possible polarizations (a, /i) = (H, V), (V, H), (H, H), {V, V), 
( + , — ), ( — , +), +), ( — , — ). To formulate the probabilities 
conditional on any source states, we only need to relate the source 
state with the beam-splitter's input state. Suppose the source state in 
photon-number space is ® Pb and 



PA=^„an\n){n\ 
PB=^„K\n){n\ 



(32) 



After some loss channel, the state changes into the beam-splitter's 
input state as Eq.(30). Suppose the transmittance for the channel 
between Alice (Bob) and UTP is 1]^ (1b)- Using the linear loss 
model of Eq. (3) we have 

/''. = E„>.. ««'7a-(1-'?J""*'C^' 

(33) 

We now arrive at our major conclusion: 

Major conclusion: Formulas of p^''(fcifc2) in the earlier subsection 
together with Eqs. (31,33) complete the model of probabilities of 
different events conditional on any source states, i.e., the gains. 
Using Eqs. (8,9), one can also model the observed error rates of 
any source states. 

E. 3-intensity decoy-state MDI-QKD. Using the Major conclusion 
above, we can model the gains and the error rates with a 3-intensity 
decoy-state MDI-QKD method^'^ ''^ We assume that Alice (Bob) has 
three intensities in their source states, denoted as 0, f.ij^, (0,^ig,/(^) . 
Denote p^ [py) as the density operator for source x {y) at Alice's 
(Bob's) side, and x {y) can take any value from 0, ji^, p'j^ (0,jig,p'g) . 



Po = |0>(0|;p,,= ^a,|/:)(fc|; p„, = ^ 4|/c)(fc|; 

k k 

P,,= Y.^k\k){k\; p,,^^Y.^[\k){k\ 



(34) 



Then we have the expression for the low bound of the yield of single- 
photon pulse pairs 



-yX ^ -yXX . 



<h'2[sl^-^)-a,b2[sl.^^,-\^ 



(35) 



a\ai(b'2bi-b2b\) 
and their upper bound of the phase flip-error rate 

pX cX _uX qX _ uX qX , pX cX 
X ^ X,U _^ti,f''^t',fi ^ftO'^ftO ^0,,i''0,,i^-^0,0''0,0 

e^i^e^i = TFx (36) 

^ 11 

With the results above, now we can calculate the key rate with the 
formula"*'«-« 



R>a'^b[Yl^[l-H 



(efi)]-S^y/(£^v)H(£fv) (37) 



Numerical simulations 

Using all the above correspondence, we can numerically simulate the 
gains and error rates of any source states. Taking as an example, we 
consider the source of a HSPS from parametric down- conversion 
processes". It originally has a Poissonian photon number distribution 
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Figure 2 | (a) The relative lower bound of 7fj and (b) the upper bound of 
efj for different photon sources. The solid lines (WO) represent the results 
of using infinite-decoy state method, and the dashed or dotted lines (Wl, 
P 1 or S 1 ) represent using three-decoy state method. Besides, W, P or S each 
corresponds to the scheme of using weak coherent sources'", possonian 
heralded single photon sources'^ or sub-possonian heralded single photon 
sources", individually. Here the superscript X represents in X basis, and at 
each point, we set /( = 0.05, and optimize the value for /('. 



when pumped by a continuous wave (CW) laser''^, written as: 



(38) 



where x is the the average intensity of the emission light. However, 
after chosen a proper gating time and triggered with a practical single 
photon detector, a sub-Poissonian distributed source state can be 
obtained, which can be expressed as: 



[P^°''d, + (l- 



X" 

JrT- 



1)! 



]|0>(0|- 



\n){n\ 



(39) 



where P*^" is the correlation rate of photon pairs, i.e., the probability 
that we can predict the existence of a heralded photon when a her- 
alding one was detected; is the dark count rate of the triggering 
detector. 

In the following numerical simulations, for simplicity, we assume 
the UTP lies in the middle of Alice and Bob, and aU triggering 
detectors (at Alice or Bob's side) have the same detection efficiency 
(75%) and the same dark count rate (10"''). We also assume all 
triggered detectors (at the UTP's side) have the same detection effi- 
ciency (they are attributed into the channel loss), and the same dark 



0.10 



N 

LU 



0.05 




0.00 

20 40 

Channel loss (dB) 

Figure 3 | (a) The gain and (b) the quantum error-bit rate of the signal 
pulses in Z basis for different photon sources. The solid lines (WO) 
represent the results of using infinite-decoy state method, and the dashed 
or dotted lines (Wl, PI or SI) represent using three-decoy state method. 
Besides, W, P or S each corresponds to the scheme of using weak coherent 
sources, possonian heralded single photon sources'' or sub-possonian 
heralded single photon sources*^, individually. Here at each point, we set 
j.1 = 0.05, and optimize the value for 
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Figure 4 j (a) The final key rate for different photon sources. The solid lines 
(WO) represent the results of using infinite-decoy state method, and the 
dashed or dotted lines (Wl, PI or SI) represent using three-decoy state 
method. Besides, W, P or S each corresponds to the scheme of using weak 
coherent sources, possonian heralded single photon sources" or sub- 
possonian heralded single photon sources*^, individually. Here at each 
point, we set ;i = 0.05, and optimize the value for 
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count rate (3X10 ''). Besides, we set reasonable value for the system 
misalignment probability Ej = 1.5%, and for the correlation rate of 
photon pairs P'^" = 0.4'=. 

Fig. 2(a) and (b) each shows the relative low bound of yfj and the 
upper bound of efj changing with channel loss for different source 
states, i.e., the weak coherent sources (W), the possonian heralded 
single photon sources (P) and the sub-possonian heralded single 
photon sources (S). The solid line represents the result of using 
infinite number of decoy state method (WO), and the dashed or 
dotted lines (Wl, PI or SI) are the results of using three-decoy state 
method. 

Similar to Fig. 2(a) and (b). Fig. 3(a) and (b) each shows corres- 
ponding values of the gains (S^/^,,) and the quantum bit-error rates 
(QBER) (-E^/^,/) of the signal pulses in Z basis for different source 
states. And Fig. 4 presents the final key rate changing with channel 
loss. 

See from Fig. 4, we find that the sub-possonian heralded single 
photon sources can generate the highest key rate at lower or mod- 
erate channel loss (s64 dB). Because within this range, its signal 
state has a lower QBER than in the weak coherent sources, and a 
higher gain than in the possonian heralded single photon sources as 
simulated in Fig. 3 (a) and (b). However, at larger channel loss 
(S64 dB), the possonian heralded single photon source shows better 
performance than the other two, this is mainly due to its much lower 
vacuum component which may play an essential role in the key 
distOlation process when suffering from lager channel loss. 

Conclusions 

In summary, we have presented a model for simulating the gains, the 
error rates and the key rates for MDI-QKDs, which can be applicable 
to the schemes of using arbitrary probabilistic mixture of different 
photon states or using any coding methods. This facilitates the per- 
formance evaluation of the MDI-QKD with phase randomized gen- 
eral sources, and thus makes it a valuable tool for devising high 
efficient QKD protocols and for studying long distance quantum 
communications. 
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